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METHOD AND SYSTEM FOR NEAR-INFRARED FLUORESCENCE 
CONTRAST-ENHANCED IMAGING WITH AREA ILLUMINATION 
AND AREA DETECTION 

TECHNICAL FIELD OF THE INVENTION 

The present invention relates generally to biomedical imaging and, more 
particularly, to a method and system for near infrared fluorescence contrast-enhanced 
imaging with area illumination and area detection. 

BACKGROUND OF THE INVENTION 

Fluorescent contrast agents have been developed and/or employed that 
enhance the optical detection of diseased tissues. These agents excite and re-emit at 
near-infrared (NTR) wavelengths, which deeply penetrate tissues. Localization of a 
contrast-enhanced target in three dimensions necessitates: (1) rapid and precise 
acquisition of fluorescence measurements; (2) accurate prediction of these 
measurements using an appropriate theoretical model of light propagation through 
tissues; and (3) tomographic reconstruction of model parameters, which include 
optical and fluorescent properties of tissues, using an optimization routine that 
minimizes the differences between the experimental and predicted fluorescence data. 
Steady-state and time-resolved experimental measurements of diffuse fluorescence 
have been used to localize a contrast-enhanced target. The steady-state measurements 
provide information regarding the spatial origin of fluorescent photons, and the time- 
resolved measurements provide additional information regarding fluorescence decay 
kinetics, which may render functional information about the environment in which the 
fluorophores reside. 

Experimental measurements of diffuse fluorescence acquired using point 
illumination and point detection techniques have been incorporated into tomographic 
reconstruction algorithms. In order to adequately probe a tissue volume using these 
techniques, illumination and detection must occur via fiber optic at several points 
about the tissue volume. Disadvantages with point illumination and point detection 
include (1) possible failure to sufficiently excite the contrast-enhanced target and (2) 
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the production of sparse data sets, which make the underdetermined inverse problem 
difficult to solve. 

SUMMARY OF THE INVENTION 

A method and system for near-infrared fluorescence contrasted-enhanced 
imaging with area illumination and area detection are provided. In one embodiment, 
a method for biomedical imaging includes directing time- varying excitation light at a 
surface area of a light scattering material, the material comprising a fluorescent target. 
Time-varying emission light from the fluorescent target is detected, substantially at a 
two-dimensional sensor surface, in response to the time-varying excitation light 
stimulating the fluorescent target. The time-varying emission light is filtered to reject 
excitation light re-emitted from the material. A three-dimensional image of the 
fluorescent target is generated based on the detection substantially at the sensor 
surface. 

Embodiments of the invention provide a number of technical advantages. 
Embodiments of the invention may include all, some, or none of these advantages. 
An approach employed herein involves area illumination and area detection on the 
same tissue-like surface. This approach overcomes the disadvantages associated with 
point illumination and point detection given that (1) fiuorophores within a large tissue 
volume are more efficiently excited via area illumination and (2) significantly more 
data is provided via area detection. In addition, non-contact imaging is possible. 
Some applications may arise in medical imaging of diseased tissues and intraoperative 
3D identification of diseased (arid tumor) lesion margin for effective and complete 
resection of diseased tissue, among others. 

Other advantages may be readily ascertainable by those skilled in the art from 
the following figures, description, and claims. 

BRIEF DESCRIPTION OF THE DRAWINGS 

For a more complete understanding of the invention, and for further features 
and advantages, reference is now made to the following description, taken in 
conjunction with the accompanying drawings, in which: 



WO 2004/113889 



PCT/US2004/019792 



3 

FIGURE 1 is a schematic illustrating a system for near-infrared ("NIR") 
frequency-domain photon migration ("FDPM") measurements according to one 
embodiment of the present invention; 

FIGURE 2 is a schematic illustrating a system for detecting re-emitted 
excitation light in accordance with one embodiment of the present invention; 

FIGURE 3 is a flow chart illustrating a method for tomographic fluorescent 
imaging in accordance with one embodiment of the present invention; 

FIGURE 4 illustrates one embodiment of generation step of FIGURE 3; and 

FIGURES 5A-C illustrate a reconstructed three-dimensional image of a 
fluorescent target. 

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS OF THE INVENTION 

FIGURE 1 is a schematic illustrating a system for excitation area illumination 
and area detection using near-infrared ("NIR") frequency-domain photon migration 
("FDPM") imaging techniques according to one embodiment of the present invention. 
In particular, system 100 generates three-dimensional ("3D") images of fluorescent 
target 102 within a light scattering material (e.g., tissue) based on emission light 
detected at a two-dimensional ("2D") sensor surface. At a high level, system 100 
delivers an intensity-modulated excitation light 104 into the tissue of a body 106 to 
cause stimulated emission of fluorescent target 102. Emission light 108 from the 
fluorescent target 102 is detected by intensified charge-coupled camera (ICCD) 112 
via gain modulated image intensifier 114. A 3D image of fluorescent target 102 is 
generated based, at least in part, on excitation light 104 and emission light 108. 
System 100 provides several advantages over the prior art such as, for example, the 
ability to reconstruct a 3D image from data collected at a 2D surface with contact. 
The present invention contemplates more, less, or different components than those 
shown in system 100 of FIGURE 1. 

Referring to FIGURE 1, system 100 includes a laser diode 116 coupled to a 
laser DC bias 118 and a reference frequency generator 120. As a result, laser diode 
116 generates light source 122, intensity-modulated excitation light. In one 
embodiment, light source 122 comprises 785-nanometer ("nm") excitation light 
provided by a 70-milliWatt ("mW") laser diode 116. Additionally, a temperature 
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controller (not illustrated) may be include in system 100 to assist in maintaining the 
optical power and the wavelength of laser diode 116; however, other suitable 
wavelengths and sources may be used.. Light source 122 may comprise other suitable 
NIR light such as, for example, 700-nrn to 900-nm light. Light source 122 passes 
through lens 124A to produce an expanded beam of excitation light 104 for 
illuminating a surface area 126 of body 106. Surface area 126, in one embodiment, is 
approximately 9 centimeters squared (cm 2 ) in size. Surface area 126 may comprise 
other suitable areas ranging from 1 cm 2 to 100 cm 2 . Surface area 126 is illuminated 
with excitation light 104 for stimulating fluorescent target 102. Fluorescent target 
102 comprises a fluorescent contrast agent (e.g., Indocyanine Green ("ICG")) that 
emits light in response to excitation light 104. In one embodiment, ICG is excited at 
780 nm and emits at 830 nm, has an extinction coefficient of 130,000 M _ l cm -1 , a 
fluorescent lifetime of 0.56 ns and a quantum efficiency of 0.016 for the 780/830 nm 
excitation/emission wavelengths. After illuminating surface area 126, light is re- 
emitted from body 106. 

Re-emitted light 110 includes emission light 108 and reflected excitation light. 
To minimize the excitation light delivered to image intensifier 114, a band-pass filter 
128 centered around the emission wavelength and a band-rejection filter 130 centered 
around the excitation wavelength are used in combination to substantially reject the 
excitation wavelength from re-emitted light 110 to provide emission light 108. 
Emission light 108 substantially comprises light emitted by fluorescent target 102 as a 
result of excitation light 104 stimulating fluorescent target 102. Band-pass filter 128 
and band-rejection filter 130 may be coupled via lens 124B (e.g., 50-mm lens). In 
one embodiment, band-pass filter 128 and band-rejection filter 130 comprise an 830- 
nm image quality bandpass interference filter and a 785-nm holographic band- 
rejection filter, respectively. Alternatively, system 100 may use any suitable 
combination of filters such as, for example band-pass filters, long-pass filters, 
holographic notch filters, band-rejection filters or any combination thereof, or any 
optical scheme to reject the excitation wavelength enabling registration of the 
emission wavelength. After the combined filters substantially remove the excitation 
wavelength, emission light 108 is amplified by image intensifier 114. 
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In one embodiment, image intensifier 114 includes a photocathode face that 
converts photons to electrons, a multi-channel plate ("MCP") that multiplies the 
electronic signal by avalanche multiplication, and a phosphorescent screen that 
converts electrons into an optical image. In one embodiment, image intensifier 1 14 is 
a fast intensifier of the variety manufactured by Litton Electronics, Inc., which 
enables modulation by applying a DC bias via image intensifier bias 132 and an RF 
signal via an RF amplifier 134 between the photocathode and the MCP; however, 
other suitable image intensifiers may be utilized. In the illustrated embodiment, the 
intensifier gain of image intensifier 114 and the intensity of excitation light 104 are 
phase-locked by a 10 MHz reference signal provided by oscillator 120. By 
modulating laser diode 116 and image intensifier 114 at the same frequency, 
homodyning results in a steady-state image on the phosphor screen. The image from 
the phosphor screen is focused through a lens 124C (e.g., 105-mm lens) onto a two- 
dimensional sensor surface 138 of camera 112, in order to provide area detection. As 
used herein, the term "two-dimensional sensor surface" means a substantially two- 
dimensional surface determined by the sensors of camera 112. Sensor surface 138 
may be flat, curved, bent, or arranged in any suitable surface. In one embodiment, 
sensor surface 138 comprises a 16 cm 2 flat surface. Camera 112, in the illustrated 
embodiment, has a 512 x 512 array of CCD detectors configured to provide a 
corresponding pixelated image. 

In operation, following each acquired image, a phase delay between image 
intensifier 114 and laser diode 116 is induced by stepping the phase of image 
intensifier 114 between 0 and 360° relative to the phase of laser diode 116 
modulation. Preferably, control between oscillator 140 and tire processor is obtained 
by a conventional general purpose interface bus ("GPIB") 142; however, other 
suitable interfaces may be utilized. Images from the phosphorescent screen image 
intensifier 1 14 may then be gathered at each phase delay by camera 1 12 to determine 
phase-sensitive fluorescence intensity. In one embodiment, each image, may contain 
128 X 128 pixels of intensity counts and be acquired by use of a 0.4-s CCD exposure 
time. A fast Fourier transform of the phase-sensitive fluorescence intensities acquired 
at each CCD pixel may be performed to compute the amplitude and phase-shift of 
excitation light 108 at each CCD pixel. Prior to generating a 3D image of fluorescent 
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target 102 based on the amplitude and phase-shift of emission light 108, the predicted 
amplitude and phase-shift of emission light 108 may be determined by the coupled 
frequency-domain diffusion equations (1) and (2) discussed below; however, the 
amplitude and phase may be determined by any other suitable mathematical model of 
light propagation. However, the solution to these coupled equations require the 
spatial distribution of amplitude and phase-shift of incident excitation light 104. 
Since it is difficult to create an excitation beam 104 of intensity-modulated excitation 
light that is spatially uniform in both amplitude and phase-shift, the spatial 
distribution of excitation light 104 may be experimentally determined. 

In one embodiment, system 100 may be altered to experimentally determine 
the amplitude and phase-shift of incident excitation light 104 in the following manner: 
(1) linear polarizers 144 are positioned at the output of laser diode 116 and the input 
of CCD camera 1 12; (2) aperture of lens 124 may be minimized; and (3) bandpass 
filter 128 and band-rejection filter 130 are replaced with a neutral density filter 146 to 
provide system 200 illustrated in FIGURE 2. Generally, re-emitted light 110 includes 
both specularly (or singly) reflected excitation light and multiply reflected excitation 
light. Since singly reflected excitation light may be representative of incident 
excitation light 104, re-emitted light 110 is compared to multiply-reflected excitation 
light to determine the spatial distribution of the singly reflected excitation light and 
thus the spatial distribution of incident excitation light 104. 

In one aspect of operation, the polarization angle of linear polarizers 144 are 
first oriented parallel to one another to predominately detect specularly reflected 
excitation light. Images of phase-sensitive excitation light intensity of this detect light 
may be determined using the procedure described above with regard to FIGURE 1, 
predominantly resulting in images representative of specularly reflected excitation 
light. In order to substantially remove multiply scattered excitation light intensities 
from these images, the polarization angles of linear polarizers 144 may be adjusted to 
determined images of phase-sensitive excitation light intensity for multiply reflected 
light. In general, linearly polarized light incident on a tissue surface and multiply 
reflected by fluorescent target 102 reemerges from the tissue depolarized with equal 
probability of having a polarization angle parallel or perpendicular to incident 
polarization angle. Therefore, to substantially remove image intensities representative 
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of multiply reflected excitation light from the images of phase-sensitive excitation 
light intensity, a second set of images of phase-sensitive excitation light intensity 
representative of multiply reflected excitation light may be acquired after orientating 
the polarization angles of polarizer 144 perpendicular to one another. The second set 
of images may be subtracted from the first set of images to provide phase-sensitive 
excitation light intensity of the specularly reflected excitation light. Afterward, a fast 
Fourier transform of the phase-sensitive images representative of the specularly 
reflected excitation light may be performed to determine the spatial distribution of the 
amplitude and phase of incident excitation light 104. Once the spatial distribution of 
both incident excitation light 104 and emission light 108 are determined, the 
mathematical expression discussed in relation to FIGURE 4 may be used to iteratively 
converge upon the optical properties of fluorescent target 102. 

FIGURE 3 is an example flow diagram illustrating a method 300 for 
generating a 3D image from light detected on a two-dimensional sensor surface. 
Method 300 is described with respect to system 100 and 200 of FIGURES 1 and 2, 
respectively, but method 300 could also be used by any other system. Moreover, 
system 100 and 200 may use any other suitable techniques for performing these tasks. 
Thus, many of the steps in this flow chart may take place simultaneously and/or in 
different orders as shown. Moreover, system 100 and 200 may use methods with 
additional steps, fewer steps, and/or different steps, so long as the methods remain 
appropriate. 

At a high level, method 300 illustrates two processes for generating a 3D 
image: an emission-light process and an incident-light determination process. The 
emission-light determination process is illustrated in steps 302 to 316, and the 
incident-light determination step is illustrated in steps 3 18 to 332. Method 300 begins 
at step 302 where a time-varying excitation light 104 is directed at a surface area 126 
of body 106 to stimulate fluorescent target 102. Next, at step 304, re-emitted light 
1 10 from body 106 is filtered by band-pass filter 128 to pass a subband of re-emitted 
light 110 including the emission wavelength to band-rejection filter 130 while 
substantially rejecting other wavelengths. Band-rejection filter 130 rejects a subband 
of the received light, including the excitation wavelength, and passes emission light 
108 to image intensifier 114 at step 306. At step 308, image intensifier 132 is phase- 
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locked to laser diode 116 and amplifies emission light 108. In order to determine 
phase-shift data for emission light 108, image intensifier 132 is stepped between 0° 
and 360° relative to the phase of laser diode modulation at step 310. Next, at step 
312, the amplified emission light 108 is detected at a two-dimensional surface area 
138 of CCD camera 1 12 to electronically recover intensity data of emission light 108. 
Images of phase-sensitive fluorescence intensities is detennined for each CCD pixel 
at step 314. At step 316, a fast Fourier transform of the phase-sensitive fluorescence 
intensities is performed to compute spatial distribution of the fluorescence intensity 
and phase-shift for each CCD pixel. In short, system 100 has now experimentally 
determined the spatial distribution of emission light 108 and is ready to determine the 
spatial distribution of incident excitation light 104. 

The excitation-light determination step begins at step 318 where linear 
polarizers 144 are installed at the output of laser diode 116 and the input of image 
intensifier 114. Next, at step 320, bandpass rejection filter 128 and band-rejection 
filter 130 are replaced with neutral density filter 146. In order to determine the spatial 
distribution for excitation light, the polarization axis of linear polarizers 144 are 
aligned along the same axis at step 322. A first set of images of phase-sensitive 
excitation intensities is determined at step 324, which predominantly represents 
specularly reflected excitation light. At step 326, linear polarizers 144 are aligned 
with their polarization axis substantially perpendicular to detect multiply scattered 
excitation light. A second set of images of phase-sensitive excitation intensities is 
determined at step 328, which predominantly represents multiply reflected excitation 
light. In order to substantially remove multiply reflected excitation intensities from 
the first set of images, the second set of images is subtracted from the first set of 
images to determine images of phase-sensitive excitation intensities for specularly 
reflected excitation light at step 330. Since specularly reflected excitation light is 
representative of incident excitation light 104, a fast Fourier transform of the images 
of phase-sensitive excitation intensities for specularly reflected excitation light 
computes the spatial distribution of excitation intensity and phase-shift for incident 
excitation light 104 at step 332. Finally, at step 334, a 3D image of fluorescent target 
102 is generated based on the spatial distribution of intensity and phase-shift for 
incident excitation light 104 and emission light 108. FIGURE 4 illustrates one 
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embodiment of determining the 3D image corresponding to step 334 above, in 
accordance with the penalty/modified barrier function (PMBF) method. 

FIGURE 4 illustrates one embodiment of step 334 of FIGURE 3 for 
generating a 3D image of fluorescence absorption coefficients. Method 400 is 
described with respect to system 100 and 200 of FIGURES 1 and 2, respectively, but 
method 400 could also be used by any other system. Moreover, system 100 and 200 
may use any other suitable techniques for performing these tasks. Thus, many of the 
steps in this flow chart may take place simultaneously and/or in different orders as 
shown. Moreover, system 100 and 200 may use methods with additional steps, fewer 
steps, and/or different steps to generate a 3D image of fluorescence absorption 
coefficients. 

To aid in understanding the various mathematical aspects of method 400, the 
following tables of selected variables is listed: 

r position 

c speed of light 

S source 

n number of nodal points 

co angular modulation frequency 

\i axi absorption coefficient due to chromophore 

\x Uxf absorption coefficient due to fluorophore 

\x a total absorption at the emission wavelength is equal to the absorption 

owing to nonfluorescent chromophores fj.^. and fluorescent 
choromophore, \i a ^ 

D optical diffusion coefficient =1/3 ^ji^ m +P-' Sx ] 

p Sjt m isotropic scattering coefficient 

[i 0x m total absorption coefficient 

£> complex fluence 

<j> quantum efficiency 



WO 2004/113889 



PCT/US2004/019792 



10 

x fluorescence lifetime 

Q volume 

I^q intensity amplitude 

9 phase-shift 

5 K global stiffness matrices 

b global vector in FEM stiffness equation containing source terms 

E error function 

N D number of detectors 

Ns munber of source 

10 g gradient vector VE 

d Newton direction 
a step length 

H Hessian matrix V 2 E 

I lower bound of the absorption coefficient owing to fluorophore 

15 n upper bound of the absorption coefficient owing to fluorophore 

X 1 Natural coordinate 

B (y, a ) active set (lower bound) 

B' (\i a ) active set (upper bound) 

C ( \i a ) free variable set 
20 s tolerance of lower and upper bound 

1 Lagrange multipliers 

X 1 Lagrange multiplies for lower bound 

X u Lagrange multipliers for upper bound 

M modified barrier penalty function 

25 ^ bound constrain 

/ [fy (^a )) quadratic penalty/barrier term of PMBF 

cp (hi (\i a )) quadratic extrapolation method 

Sj shifting parameter 
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barrier parameter 



P 



parameter defines how close the extrapolation will get to singularities 
of logarithm term 




coefficient of the quadratic function 



y 



reduction factor of the barrier parameter 



* 



outer iteration number 



J 



inner iteration number 



yf ' 72 » d\ and d 2 convergence criteria 



reconstruction error 



Generally, method 400 uses the modified barrier penalty function to 
reconstruct a 3D image of fluorescence absorption coefficients. However, prior to 
discussing method 400 additional information will be provided. 

Diffusely propagated fluorescence is predicted by the coupled diffusion 
equations (1) and (2) given the optical properties of the tissue. The numerical 
solution of the diffusion equation may be needed to provide sufficiently accurate 
estimates of the fluence on the boundary of the tissue for the image reconstruction (or 
inverse) problem. The coupled diffusion equation used in optical tomography at a 
given modulation frequency /(co = 2nf rad) is given by: 



-V{D X (?) V<D, (r.co)] + jjp+u % . (?) (?) k (?,a>) = S 
-V-[D m (r)VO m (f ) a>)]+ l * +M (?)|<D W (?,*>)= # -±— < 



on Q (1) 



(2) 



Here, fl^and 0,„ are the AC components of the excitation and emission photon 
fluence (photons/cm 2 /sec) at position r, respectively; a is the angular modulation 
frequency (rad/s); c is the speed of light within the medium (cm/s); u^. is the 
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absorption due to chromophores (cm -1 ); and n 0j/ is the absorption due to 
fluorophores (cm' 1 ). D x>m is the optical diffusion coefficient equivalent to 
\j4^i axm 4- Vs xm )> wnere Vsxn i s ^ e isotropic scattering coefficient (cm" 1 ) and 
Ha^is the total absorption coefficient, at the respective wavelengths. The total 
absorption at the emission wavelength, p. flm , is equal to the absorption owing to 
nonfluorescent chromophores, n % . , and fluorescent chromophores, ]x„ mf . The right 

hand term of equation (2) describes the generation of fluorescence within the medium. 
The term ij> represents the quantum efficiency of the fluorescence process, which is 
defined as the probability that an excited fluorophore will decay radiatively, and t is 
the fluorophore lifetime (ns). Note that the source term couples with the solution of 
excitation fluence described by equation (1). The numerical solutions for the 
excitation and emission fluence distributions given in equations (1) and (2), 
respectively, may be obtained using the Robin boundary condition. Other suitable 
mathematical models of light propagation and emission generation may be used to 
obtain the excitation and/or emission fluence distributions. 

The source term iS in equation (1) depends on the measurement types. For 
point illumination, the excitation source is represented as 

S= n is5(r-r s ) (3) 

where r s is the positional vector of the illuminating point on the surface of the 
phantom; 5 is the Dirac delta function representing the source of excitation at a single 
point; «j is the total number of simultaneous points of illumination, which is typically 
one. For area illumination, the area of the planar source can be discretized into a 
number of triangular finite elements and the source within a surface element in area 
coordinate is represented by 

s=iASi (4) 
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where the!,] are the natural co-ordinates for the triangle and S { is the source value at 
each node (i). The source term corresponding to the area illumination'' can be 
approximated as the summation of simultaneous and closely position multiple point 
sources as given by equation (4). 

The Galerlcin finite element method with tetrahedral in three-dimensions may 
be used to solve the differential equations (1) and (2). Local stiffness matrices may be 
formed for each element and then these matrices may be combined into two global 
stiffness matrices, K x for the excitation equation and K m for the emission equation. 
The solutions of the differential equations are obtained by solving the complex 
equations: 

K ,, m = (5) 

where b xm is the right hand side of excitation and emission equations. The 
discretization process converts the continuous problem into a problem with finite 
number of unknowns by expressing the unknown field variables, which are the 
absorption coefficients owing to fluorophore \i a ^ in terms linear interpolation 
functions within each element. The linear functions are defined in terms of the values 
of the field variables at each node of the element. Therefore, the nodal values of the 
field variables become the unknowns to the image problem. 

From the solution of <S> x m , the values of the amplitude, Jac xh . and phase 
shift, Q x m , can be determined from the relationship: 

® X ,m=lAC XJtt exp{-®x,m) (6) 

Two finite elements meshes may be employed when solving the modified 
barrier penalty function. The first may contain 33431 nodes with 163285 tetrahedral 
elements (denoted below as Mesh 1) while the second, more refined mesh may 
contain 38919 nodes with 193137 tetrahedral elements (defined as Mesh 2). 1205 
seconds and 1738 seconds of CPU time may be required to solve the coupled 
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excitation and emission equations for Mesh 1 and Mesh 2, respectively. All the 
computations in this paper may be performed on a SUN-blade 100 UNIX workstation. 

Information regarding the error function will also be provide prior to 
discussing method 400. The error function considered for the image reconstruction 
problems is as followed (for brevity we will use \i a for henceforth): 

E(^h\ N i^log[z p ) c - log(z p )^ [log{z p l - log(z p f m ^ (7) 



where E is the error function associated with measurements, Z p is comprised of 
referenced fluorescent amplitude, lAC re j > ^ referenced phase shift 9,.^ 

measured at boundary point, p, in response to planewave illumination. Specifically 
the referenced measurement at boundary point p is given by: 

> Zp=lAC refp exp(i%ef p ) (8) 

The experimental values of fluorescence modulated amplitude J ACp , may be 
divided by the maximum fluorescence experimental value from the collected 
reflectance data set, lAC max ^ n order to compute the referenced amplitude values, 
lAC ref The experimental fluorescence phase shift, at the position where the 
maximum fluorescence lAC max was obtained, may be determined as B at j 4r and 
subtracted from all experimental fluorescence phase shift values, 0^ in the collected 
reflectance data sets to obtained referenced phase shift values, ® re f p ■ These represent 
the "referenced" measurements: 
I AC 

^max 
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The predicted fluorescence and Q refp may be computed by solving the 

coupled diffusion equations (1 and 2) using the finite element method and may be 
similarly normalized to obtain the calculated values. Other suitable methods may be 
used to calculate these values. The referenced data are normalized in this manner in 
order to reduce, minimize, or eliminate instrument functions. 

The error in equation (7) is reported for values of tissue absorption owing to 
fluorescent contrast agent. The subscript c denotes the values calculated by the 
forward solution and the subscript m denotes the experimental value. N D denotes the 
number of detectors on the surface. The superscript * denotes the complex conjugate 
of the complex function of the referenced measurements. 

The discussion regarding the coupled diffusion equations (1) and (2) and the 
error function (7) provides a basis for the modified barrier penalty function, M 
(termed hereafter as the barrier function), which is stated as followed: 



(11) 



where r\ is a penalty/barrier term, n is the number of nodal points. X 1 and X u are 
vectors of Lagrange multipliers for the bound constraints for the lower and the upper 
bounds, respectively. From equation (11) one can see that the simple bounds are 
included in the barrier function M. 

Function,/ may be defined as: 

f(^ a )h\ l0g [ Si for r*<W * (12) 

MkM) fori=l,2.:.,nifh i ( } x a )<-ps i r\ 

where 

<tfaba)) = y<kM? +ff?ftrOi fl ))+ft 3 (13) 
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hi{\i a ) represents all the bound constraints {(u a; - /) and(« - \i ai )}; and s t is 
positive-valued shifts used in the enlargement of the feasibility region for the lower 
and the upper bounds, respectively. In one embodiment, the s t is set to 1.0 for all 

computations. The coefficients qj, qf and qf of the quadratic function in equation 
(13) are chosen so that their values, and their first and second derivatives match those 
of the logarithmic term in equation (12) at the point ^(u*) = -pj,q. The coefficients of 
the quadratic as given as: 



1i = 



kn(i-P)) 2 

l-2p 

3(2 



(14) 
(15) 



[2(l-p) : 



pj + fcfoG-pj) (16) 



The coefficients are dependent on both p and T] , and thus have to be updated for any 
changes of p and t] . The values of P defines how close to the singularities the 
logarithmic terms described in equation (12) are extrapolated. Consequently we 
define a very simple differentiable modified barrier penalty function by equation (1 1) 
that is used to assess the quality of the solution. This function has an extremely 
simple structure so that the overhead for using the barrier function instead of the 
original error function given by equation (7) is negligible for function evaluation or 
for computing gradients and Hessian as discussed below. 

The PMBF method 400 is performed in two stages in, an inner and an outer 
stage illustrated in FIGURE 4. The outer iteration updates the Lagrange multiplier X, 
parameter p, and the barrier parameter T|. Formulations used to update these 
parameters at each outer iteration are provided in the discussion below. Using the 
calculated values of the parameter X, an unconstrained technique may be applied to 
minimize the penalty barrier function described by equation (1 1). Once satisfactory 
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convergence within the inner iteration is reached, the variables in the outer iteration 
are recalculated to check for satisfactory convergence in the outer iteration. If 
unsatisfactory, the Lagrange multipliers, the barrier parameters are updated and 
another unconstrained optimization is started within the inner iteration. 

Turning to method 400, method 400 begins at step 402 where initial values of 
H a , r\ and p are selected. The initial value for u fl , the absorption coefficient due to 
fluorophore, is a positive value typically selected from the range of 0 to 0.1 cm" 1 . The 
initial value for n, the barrier parameter, is an arbitrary positive value and typically 
10" 1 or 10" 2 . The initial value for (3 is a parameter that defines how close to the 
singularities the logarithmic terms in the modified barrier penalty function are 
extrapolated. The values of p are typically not too small (e.g., 0< p<0.5) or very close 
to one (i.e., p = 0.99). For example, the value of p may be set to 0.9. Next, at 
decisional step 404, if the initial values of X 1 , an initial Lagrange multiplier for a 
lower bound, and X u , an initial Lagrange multiplier for an upper bound, are selected, 
then, at step 406, the values are set to equal values such as, for example, 1, 10, 100 
and 1000. If the initial values of X 1 and X u are not initially selected, then, at step 
408, the values of X 1 and X u are determined. In one embodiment, the initial values 
of X 1 and X" are determined using a vector X° that minimizes the least-squares 
unconstrained problem: 




(17) 



this gives the formula: 




(18) 



where 



^)=(v„„/'(*k)K./(4:))r 



(19) 
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The formulation of equation (18) requires calculation of all constraint gradients and 



To reduce the computational cost, only the diagonal terms of the equation (18) maybe 
used. Hence there is no need to invert the matrix. Using this least-square procedure 

the initial Lagrange multipliers for image reconstruction problem may be computed. 

After the initial values of u, a , n, (3, X 1 and X u are set, the indices k, which identifies 

the current step of the outer stage, is set to zero at step 410. 

Method 400 now executes the inner stage of the PMBF method by setting the 

indices j, which represents the current iteration of the inner stage, to zero at step 412. 

Next, at step 414, the barrier function, the gradient of the barrier function and the 

Hessian of the barrier function are calculated. The barrier function may be calculated 

by equation (1 1). During this calculation, an unconstrained optimization method may 
be used to optimize the unconstrained barrier function of equation (11) for each set of 
parameters X k and r\ k . The truncated Newton method with trust region (TN) may be 
used for the unconstrained optimization problem. The gradient based truncated 
Newton method is efficient and suitable for large-scale problem. The gradient of the 
barrier function may given by 



The gradients of the error function V^E may be calculated by the analytical method 

using reverse differentiation technique. The Hessian of the barrier function may be 
determined as followed: 



the gradients of the error function E\\i° a jand inversion of the matrix ^(nj) • The 
inversion of this large (n x n) matrix in the reconstruction problem can be very costly. 
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i?(u. a ),d maybe calculated by the finite difference scheme, 

Hba)-d=\WPa+8d)]-gfoa)} (22) 

where £(u a ) = g(p, fl ) is the gradient; = #fiu fl ) is the Hessian of the 

error function £(n a ); is the Newton direction; and 6 is (approximately) the 
square root of the machine precision. It may not safe to perform finite differencing 
directly to V 2 M because of the singularity of the logarithmic function, the final 



the formulas above . 

If the absolute value of the gradient of the barrier function is less than 10" 5 at 
decisional step 416, then the inner stage ends and execution proceeds to step 430 of 
the outer stage. Otherwise, V^M^)^' = V^ilf (fijj) is solved at step 418 to 

determine the Newton direction d. The Newton direction d is computed as an 
approximation solution of the Newton equations: 

(V 2 M)? = -VM (23) 
where V 2 M = V 2 M([i a ) is the Hessian matrix of second derivatives at the current 
solution designated by a vector \i a . The approximate solution of equation (23) may 
be obtained by the linear conjugate-gradient method but the solution may alternatively 
be determined by any suitable method. This iterative method is "truncated" before the 
exact solution is obtained. The conjugate-gradient method (CG) corresponds to 
minimizing the quadratic model 
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minQ{d) = -d T l 2 Md + d T VM 



(24) 



s.t. I ~H a ^d<u-\i n 



(25) 



as a function of d over a sequence of subspaees of increasing dimension. Criteria must 
be determined and applied to find out whether convergence has been achieved. This 
method has two levels of iterations, with the 'outer-level' iteration determining the 
step length, and the 'inner-level' iteration computing the Newton direction using 
conjugate gradient method. A further issue for unconstrained minimization is what to 
do when the current Hessian V 2 M is not positive definite. Several approaches are 
used, among the most popular of which are so-called modified Newton method in 
which dj satisfies the linear system, 



where R is a positive semidefinite matrix. Different techniques have used to calculate 
R so that (v 2 m(;c A )+ i?) is positive definite. An alternative globalization strategy in 
unconstrained optimization embodied in trust region methods, is to choose the step 
d J to minimize the local Taylor-series model subject to a limitation on the size of 



difficulty, the conjugate gradient method with trust region may be used for 
unconstrained optimization problems. The procedure is described in the following 
paragraph as followed: 



the Newton direction and the conjugate direction y 1 , The update scheme is 




(26) 



\d\ , where various norms may appear in the trust-region constraints. To alleviate this 



CG methods generate sequences p l , y l 



, where p l approximates iteratively 



P 



i+l 



= P i + aV 



(27) 



y i+1 =r i+1 



+ P y 



i 



(28) 
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where A' minimizes the quadratic function Q(d) along the direction y i , and p' is 
chosen so as to maintain conjugacy between the directions y l , The CG tenninates 
when the residual r l = - VQi[p l ) is in norm below a prescribed tolerance or when a 
negative curvature direction tj =[s 1 ^ V 2 Mj z <, 0 is encountered. This algorithm 
does not stop, but computes > 0 so that |p z + Tj/J = A, where A is the trust 
region. Then update rule is d l+1 = d l + x x y l . The algorithm stops if the residual 
r'J| is below a given tolerance, or if (y'f^My 1 = 0 . 

Next, at step 420, a new estimate of the absorption owing to fluorophore is 
detennined by [i J a +1 = \l{ + ad j , where a is the step length. The step length may be 
determined by an Armijo line search algorithm as disclosed in M. S. Bazaraa, H. D. 
Sherali, and C. M. Shetty, Nonlinear programming theory and algorithms, John 
Wiley & Sons Inc, New York., 1993. The Newton direction d may be required to 
satisfy d T VM < 0 (i.e., it is a descent direction for M at the point \i a ). The step 
length a > 0 may be chosen to guarantee that M(p. a ) < M(jj. a ), along with Wolfe's 
conditions disclosed in P. Wolfe, "Convergence condition for ascent method," SIAM 
Rev., 11, 226-253, 1969, which are designed to guarantee convergence to a local 
minimum. 

At decisional step 422, its is determined whether is an active variable. 
In general, if the variables of an unconstrained problem are very close to a bound in 
the inner iteration, then these variables are fixed (call active variables). Hence the 
start of the next unconstrained optimization uses the free variables only. The 
variables in B(u. a ) and B'()a a )are called active variables and the variables in C(jj. a ) 
are called free variables 



B(p, a ) = {i: l<\x at <l +s} 
B'dx a )={i:u -s<u. fl . <u } 



(25) 
(26) 
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C(fi fl )={/ •■/ + s^n fl/ <« -e } (27) 

where z is the number of nodal points i=l,...n , I is the lower bound of the absorption 
coefficient owing to fluorophore, « is the upper bound of the absorption coefficient 
owing to fluorophore, and the tolerance s should be sufficiently small, such that 

/ + s < u - s (28) 

For example, the tolerance s may be set to 10" 3 , 10" 4 , 10" 6 , or other suitable values. 
All variables in the active set are fixed during the iteration, while the other variables 
in set Care free. After each iteration, the active variables may be updated from 
current free variables, i.e., variables may be found that lie within the bounds and add 
to previous active variables. The active sets B and B' are updated at the inner 
iteration and the process is continued until the gradient of the barrier function M 
described by (1 1) can not be reduced further with the current free variables. 

The variables \i a . in C= jn a . | / < n«. < w} have three mutually exclusive 



(i) 0<v a .-l<V llai M(yi a Xr l ) (29) 

(ii) V^.M(n fl ,M)<>i fl .-«<0 < 3 °) 
(hi) ^-^V^.M^Ati)^^.-/ (31) 

The vector of u fl . that satisfies (i) or (ii) are now the dominated variables. The ones 
satisfying (i) are said to be dominated above; and those satisfying (ii) dominated 
below; and those satisfying (iii) are the floating variables. 

If a convergence condition is satisfied at decisional step 424, then inner stage 
ends and method 400 proceeds to step 430 of the outer stage. In the illustrated 
embodiment, the inner iteration can be considered fully converged if either of the 
following two conditions are satisfied. The first convergence condition is 
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Jk 



(32) 



where k is the number of the subproblem and j is the number of inner iteration. 



continued. The second convergence condition is satisfied if the fractional change <f 2 
given by (33) of the barrier function M of equation (1 1) of consecutive inner iterations 
is less than zq = 10 -5 then the inner iteration is stopped and the outer iteration are 
continued. 



If neither convergence condition is satisfied at decisional step 424, then, at step 426, 
the inner iteration indices / is set to If the inner iteration indices j does not 
satisfies the specified number of inner iterations (e.g., 5, 10, etc) at decisional step 
428, the inner iteration returns to step 418. Otherwise, the inner iteration ends, and 
execution returns to the outer iteration at step 430. 

At step 430, the outer iteration indices k is set to k+1. If the outer iteration 
indices k satisfies the specified number of outer iterations (e.g., 6) at decisional step 
432, then, at step 434, the absorption coefficient due to fluorophores \i a is set to the 

current estimated value and execution ends. If the outer iteration indices k does 
not satisfy tire specified number of outer iterations at decisional step 432, then 
execution proceeds to decision step 436. If a convergence criteria is' satisfied at 
decisional step 436, then execution ends. The outer iteration of the PMBF can be 
considered fully converged if it satisfies any of the following conditions: 



If x\ is less than 10 5 then the inner iteration is stopped and the outer iteration are 




(33) 
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71 = ^ i.o + ^| 



1.0+||iJ +1 | 



(34) 



(35) 



where k is the number of subproblem. If q is less than 10" 5 , then the inner iteration is 
stopped. If the fractional change d x given by equation (36) of the barrier function M 
given by equation (11) for the consecutive outer iterations is less than s 0 = 10~ 5 , 
then the outer iteration is stopped. 

If the outer iteration convergence criteria are not satisfied at decisional step 436, then, 
at step 438, X 1 , X u , and p are updated. The Lagrangian multipliers may be updated 



1 r\*Si + kM 



r\ k 4{q}h i ( i i a )+ qf) for i = \,2,...,n ^(fij <-pVl* 

(37) 



where X t represents all the Lagrange multipliers for the bound constraints (X 1 and 
X u ) and ^ (ix a ) represents all the bound constraints, q may be updated by a simple 
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where y > 1 is defined by the user, and normally y = 2 or y = 10, and k is the 
number of the outer iteration. 

Next, at step 440, n a . is updated at nodal point i, after k th outer iteration as 
follows: 



m = 



l if 0 < (fil- l \-l< (V^M*- 1 ),. 
u if (v M M k ~ l \z (tf* ).-/< 0 (39) 
[(x* -1 ). otherwise 



Execution returns to step 412 to return to the inner stage. 

FIGURES 5A-C illustrate x-z views of a 3D distribution of a fluorescent target 
which provides example data generated according to the teachings of the invention. 
These reconstructions used the following initial values =0.03 cm' 1 and 

A°=1000, ifij =0.003 cm' 1 and 1°=1000, and $ =0.003 cm' 1 and A,°=1000, 
respectively. A tissue phantom consisting of an 8 x 8 x 8 cm 3 cube filled with 1% 
lipid solution to mimic the scattering properties of tissue was used in the system of 
FIGURE 1 and 2. The target to be reconstructed in these initial studies was a 1 cm 3 
cube filled with 1% lipid and 1 yM Indocyanine Green solution placed 1 cm from the 
imaging surface. The top surface of the tissue phantom was illuminated with an 
expanded beam of 785-nm excitation light which was intensity modulated at 100 
MHz. The surface area illuminated was approximately 9 cm 2 . The emitted 
fluorescent signal in this area was imaged using a gain modulated intensified charge 
coupled device (ICCD) operated under homodyne conditions to enable determination 
of the amplitude and phase-delay of the light emitted on the top of the surface. Using 
a combination crossed and parallel polarizers on the incident and collected light, the 
spatial distribution of the amplitude and phase delay of the incident excitation light 
was determined. Next upon placing an 830 nm interference and imaging holographic 
filters at the photocathode of the intensifier, the amplitude and phase delay of the re- 
emitted fluorescent light owing to the submerged target was determined across the 
surface of the phantom. 
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One advantage of this procedure is that optical fihers are not required to 
illuminate or collect the light signal unlike the traditional point illumination/collection 
measurement schemes used in optical tomography. An additional advantage of this 
procedure is the ability to recontruct 3D images from data collected at a 2D sensor 
surface. 

Although the present invention has been described in detail, various changes 
and modifications may be suggested to one skilled in the art. It is intended that the 
present invention encompass such changes and modifications as falling within the 
scope of the appended claims. 
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WHAT IS CLAIMED IS: 

1 . A method of imaging, comprising: 

directing time-varying excitation light at a surface area of a tissue, the tissue 
comprising a fluorescent target; 

modulating the intensity of the excitation light at a radio frequency; 

detecting, substantially at a two-dimensional sensor surface, time-varying 
emission light from the fluorescent target in response to the time-varying excitation 
light stimulating the fluorescent target; 

filtering the time-varying emission light to reject excitation light re-emitted 
from the material; 

determining a spatial distribution of amplitude and phase-delay of the detected 
time- varying emission light; 

detecting time-varying excitation light specularly reflected by the fluorescent 

target; 

determining a spatial distribution of amplitude and phase-delay of the time- 
varying excitation light incident on the surface area based on the reflected time- 
varying excitation light; and 

generating a three-dimensional image of the fluorescent target based on the 
spatial distributions of the incident time- varying excitation light and the time- varying 
emission light. 
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2. A method of imaging, comprising: 

directing time-varying excitation light at a surface area of a light scattering 
material, the material comprising a fluorescent target; 

detecting, substantially at a two-dimensional sensor surface, time-varying 
emission light from the fluorescent target in response to the time-varying excitation 
light stimulating the fluorescent target; 

filtering the time-varying emission light to reject excitation light re-emitted 
from the material; and 

generating a three-dimensional image of the fluorescent target based on the 
detection substantially at the sensor surface. 

3. The method of Claim 2, wherein the fluorescent agent comprises 
indocyanine green ("ICG"). 

4. The method of Claim 2, wherein the surface area is greater thai 1 
millimeter squared (mm 2 ). 

5. The method of Claim 2, further comprising determining a spatial 
distribution of amplitude and phase-delay of the detected time-varying emission light. 

6. The method of Claim 5, further comprising: 

detecting time-varying excitation light specularly reflected by the fluorescent 
target; and 

determining a spatial distribution of amplitude and phase-delay of the time- 
varying excitation light incident on the surface area based on the reflected time- 
varying excitation Light. 

7. The method of Claim 6, wherein in the three dimensional image of the 
fluorescent target is based on the spatial distributions of the incident time-varying 
excitation light and the time- varying emission light. 
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8. The method of Claim 2, wherein the time-varying excitation light 
comprises a planar wave. 

9. The method of Claim 2, wherein the material comprises tissue. 

10. The method of Claim 2, further comprising modulating the intensity of 
the excitation light at a radio frequency. 

11. The method of Claim 6, wherein detecting time-varying excitation 
light singularly reflected by the fluorescent target comprises: 

detecting excitation light emitted from the light scattering material; 
detecting multiply scattered excitation light emitted from the light scattering 
material; and 

subtracting the detected multiply scattered excitation light from the detected 
excitation light emitted from the light scattering material to determine the specularly 
reflected excitation light. 

12. The method of Claim 2, wherein the two-dimensional sensor surface is 
approximately 16 cm 2 . 
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13. A system for imaging, comprising: 

a laser diode operable to direct time-varying excitation light at a surface area 
of a light scattering material, the material comprising a fluorescent target; 

an image intensifier operable to detecting, substantially at a two-dimensional 
sensor surface, time-varying emission light from the fluorescent target in response to 
the time- varying excitation light stimulating the fluorescent target; 

one or more optical filters operable to filter the time- varying emission light to 
reject excitation light re-emitted from the material; and 

an imaging apparatus operable to generate a three-dimensional image of the 
fluorescent target based on the detection substantially at the sensor surface. 

14. The system of Claim 13, wherein the fluorescent agent comprises 
indocyanine green ("ICG"). 

15. The system of Claim 13, wherein the surface area is approximately 9 
centimeters squared (cm 2 ). 

16. The system of Claim 13, the imaging apparatus further operable to 
determine a spatial distribution of amplitude and phase-delay of the detected time- 
varying emission light. 

17. The system of Claim 13, further comprising a plurality of linear 
polarizers operable to pass reflected excitation light in a first configuration and 
multiply-reflected excitation light in a second configuration; and the image apparatus 
further operable to determine a spatial distribution of amplitude and phase-delay of 
the time-varying excitation light incident on the surface area based on the reflected 
and multiply reflected excitation light. 

1 8. The system of Claim 1 7, wherein in the three dimensional image of the 
fluorescent target is based on the spatial distributions of the incident time-varying 
excitation light and the time-varying emission light. 
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19. The system of Claim 13, wherein the time-varying excitation light 
comprises a planar wave. 

20. The system of Claim 1 3, wherein the material comprises tissue. 

21. . The system of Claim 13, wherein an intensity of the time- varying 
excitation light is modulated at a radio frequency. 
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FIG. 3 2/4 ^300 
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